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Abstract 

We carried out three-dimensional simulations, with about 1.4 x 10 6 par- 
ticles, of phase segregation in a low density binary fluid mixture, described 
mesoscopically by energy and momentum conserving Boltzmann-Vlasov equa- 
tions. Using a combination of Direct Simulation Monte Carlo(DSMC) for the 
short range collisions and a version of Particle-In-Cell(PIC) evolution for the 
smooth long range interaction, we found dynamical scaling after the ratio of 
the interface thickness (whose shape is described approximately by a hyper- 
bolic tangent profile) to the domain size is less than ~ 0.1. The scaling length 
R(t) grows at late times like t a , with a = 1 for critical quenches and a = ^ 
for off-critical ones. We also measured the variation of temperature, total 

particle density and hydrodynamic velocity during the segregation process. 
PACS numbers: 64.75. +g, 68.10.-m, 47.70.Nd 
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The process of phase segregation through which a system evolves towards equilibrium 
following a temperature quench from a homogeneous phase into a two phase region of its 
phase diagram has been of continuing interest during the last decades |]J , but many problems 
still remain to be solved. This is particularly so for ffuids,when particle, momentum and 
energy densities are conserved locally; these are currently the focus of both numerical studies 
and micro-gravity experiments [§,[|. 

In this paper we present computer simulations of spinodal decomposition in a three- 
dimensional mixture of two kinds of particles that we label 1 and 2 using a novel microscopic 
dynamics and computational scheme. The particles interact with each other through short 
range interactions modeled here by hard spheres having the same mass m and diameter 
d. Particles of different kinds interact also through a long range repulsive Kac potential, 
V(r) = 7 3 £/(7r). The equilibrium properties of such a system are well understood, there is 



even a rigorous proof of a phase transition at low temperatures to an immiscible state |T0 
which in the limit 7^0 |Tl[, is described by mean field theory. When the density n is low 
enough, nd 3 -C 1, and the potential sufficiently long ranged, n^ 3 ^> 1, the free energy of 
the system is well approximated by F = ksT J[ni(f) lnni(r) + ri2(f) lnn2(r)]dr + J V(\fi — 
^2 1) ^1(^1)^2(^2)^1*^2 and the 7 — > critical temperature, which should be an upper bound 
for T c 7 at 7 > 0, is given by k B T° = \n$ U(r)df. In this regime the dynamical evolution of 
the system should be well described by two coupled Boltzmann - Vlasov equations: 

^ + ^-§4 + --§4 = J[/,/i + /2], 2=1,2 (1) 
ot or m Ov 

where fi(r,v,t) are the one-particle distribution functions, Fi(r,t) = — V JV(\f — 
r/\)nj(r/)dr/, rij{rl) = J fj(f/,v,t)dv with i,j = 1,2,2 7^ j, and J[f, g] is the Boltzmann 
collision operator for hard core interactions WM. Kinetic equations of this type have been 



proposed in |fL2|] , and if the system is quenched inside the coexistence region they will de- 
scribe gas-gas segregation |13[ into two phases, one rich in particles of type 1 and the other 
rich in particles of type 2. (Examples of gas mixtures that have a miscibility gap are helium- 



hydrogen, helium-nitrogen, neon-xenon etc. |pL3| .) We believe that the model contains the 



essential features of phase separation in general binary fluid mixtures. 

To simulate our system we modeled the Boltzmann collisional part using a stochastic 
algorithm due to Bird known as Direct Simulation Monte Carlo (DSMC), while for the 
Vlasov part we used the particle-to-grid-weighting method, well known in plasma physics 
In the DSMC method the physical space is divided into cells containing typically tens 



of particles. The main ingredients of this procedure are the alternation of free flow over 
a time interval At and representative collisions among pairs of particles sharing the same 
cell. In the particle-to-grid-weighting algorithm the particle densities are computed on a 
spatial grid through some weighting depending on the particle position, then the Vlasov 
forces are calculated on the same grid. Finally, the forces at the position of each particle 
are interpolated from the forces on the grid. The coupling of these methods, which have 
been extensively used individually, made possible our simulations of phase segregation with 
1.4 x 10 6 particles, with only modest computational resources: a typical run took about 
32 CPU hours on a 233 MHz Alpha Station. It appears that this method can be extended 
to the study of the effects of phase segregation on inhomogeneous hydrodynamical flows of 
practical importance [ITT . 



Since one of our main interests was the late time hydrodynamical regime, a delicate 
balance had to be struck between the size of the system, the range of the potential, the 
temperature and the particle density, making sure that each of the methods is used within 
its range of validity and that their combination remains computer manageable. On the 
one hand the potential must be reasonably long ranged so that the Vlasov description is 
physically appropriate and numerically sound, and on the other hand it must have a range 
much smaller than the size of the system. This restriction made necessary the use of two 
spatial grids: a somewhat coarse one for the collisions and a finer grid for the long range 
potential. It also imposed the use of quadratic spline interpolation for the calculation of 
grid quantities and a ten-point difference scheme for the calculation of the forces [[RJ . 

Our results were obtained using a system with 1382400 particles in a cube with periodic 
boundary conditions. We also studied smaller systems to identify unavoidable finite-size 
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effects. The interaction potential used was gaussian, U(x) = air~2e~ x , a > 0, but there 
is no reason to believe that different repulsive potentials would qualitatively change the 
results. All quenches were performed at a total particle density nd 3 ~ 0.01 and an initial 
temperature To, To/T c ° = 0.5. The initial conditions for each run were random positions for 
all particles and velocities distributed according to a maxwellian of constant temperature. ( 
In the DSMC evolution, as in the Boltzmann-Vlasov equations, the hard cores only enter 
in determining the collision cross sections.) The total energy of the system was very well 
conserved by the dynamics. This meant that the kinetic energy and hence the temperature 
increased as the system segregated, but at late times it changed very slowly on the time 
scale of our simulations. We indicate the final temperature T in the figures. The effective 
number of particles in the range of the potential was about 100-500. 

In the following we compare results of our simulations with available theoretical and 
experimental work and check various assumptions made in the former, e.g. the neglect of 
density and temperature variations. We are also currently investigating both formal and 
rigorous Chapman-Enskog and Hilbert expansion methods for derivation of macroscopic 



evolution equations for this model [18| . The units used for lengths and times are the mean- 
free path, A = (2^nnd 2 )~ 1 , and mean-free time, r = A/c, where c = (2/c B T /m)5 and 
T = |T C ° is the initial temperature. 

We first present results for critical quenches (equal volume fractions of the two species). 
Three different potential ranges were used, and for each one of them 10-12 quenches were 
performed. The domain size was probed using the pair correlation function, C(r,t) = 
V -1 < + f,t) >, where = (n\ — n2)/{n\ + n 2 ) is the local order parameter, 

Tlx and n 2 are the local particle densities, V is the volume and the average is over the 
different runs. We determined C(f,t) by an inverse Fourier transform of the structure 
function S(k,t) = \(f)(k,t)\ 2 ; the order parameter <f)(r,t) and its Fourier transform <p(k,t) 
were computed on a 64 x 64 x 64 cubic grid. The first zero of the spherically averaged 
correlation function, C(r,t), was used as a measure of the typical domain size, R(t). The 
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data are averages of the independent runs. 

Assuming the existence of a single characteristic length scale, the dynamical scaling pre- 
diction for the late time spherically averaged correlation function is C(r,t) ~ C(r/R(t)). 
In Fig. [I] the correlation function for potential range 7 _1 = 0.4 is plotted starting at t = 160, 
showing that the system is well within the scaling regime. 

Simple dimensional analysis of the hydrodynamical evolution equations in the limit of 
large domain sizes, appropriate for late times , yields a linear growth law for the domain 
sizes [jTQfl , R(t) oc (a/r])t, when R is below Rh = V 2 /P a an d a t% law, R(t) oc (a/p)H%, 
for R above Rh [f2^ , p3|l ; a is the surface tension coefficient, r\ is the shear viscosity and p 
is the density. We measured a directly using Laplace's law and r\ by studying the decay of 
a sinusoidal velocity profile. The time evolution of R(t) in our simulations is at late times 
(see Fig. 2) R(t) = a + b(a/rj)t, with b ~ 0.13 ± 0.2(cx = 260 and rj = 3250 for 7" 1 = 0.4, 
T/T c ° = 0.6). This numerical factor is similar to the one observed in experiments || and 
recent large scale molecular dynamics simulations 0. Whether or not a crossover to a £5 
growth occurs at later times/larger domain sizes cannot be decided by our present results. 
We estimated that we would need a system with at least 4 times as many particles as the 
present one to be able to observe the growth of domains with sizes bigger than R^. 

The linear regime starts around the time when dynamical scaling begins to hold, i.e. at 
a domain size of about 12-15 times 7 _1 , the range of the potential. This is in agreement 



with Siggia |L9[ who argued that the linear regime is due to surface tension driven flows, 
so the interfaces should be well defined, i.e. their width should be small compared to the 
domain size. This width is usually taken to be of order £, the correlation length of order 
parameter fluctuations in the bulk phases. One expects £ to be roughly equal to 7 -1 , far 
away from the critical temperature, as we are |20| . 



We also looked directly at the wall profiles separating different phases and found that 
they are approximately described by solitonic solutions |^T[. These are particle densities 
rii{z) depending on a single spatial coordinate, such that fi(z,v) = ni(z)exp(—mv 2 /2k B T) 
are stationary solutions of Eq.(l) and ni(±oo) = 7^2(^00) are the mean-field equilibrium 



densities at temperature T. The order parameter profiles which satisfy the equation and 
the ones observed in the simulations(see Fig. |j) have both approximately the hyperbolic 
tangent form pofl , tanh(z/2£ i ), with z the coordinate perpendicular to the domain wall and 
£ a parameter which characterizes the interface thickness. We found that this thickness 
is about 50 percents bigger than 7 _1 . The total density is about 20 percents smaller (for 
7 _1 = 0.4, T/T c ° = 0.6) at the interface than in the bulk, in very good agreement with 
the solitonic solution (see Fig. ^). The solitonic profiles were calculated at an effective 
temperature T e ff for which the asymptotic values of the order parameter matched the ones 
observed in the simulations. These asymptotically matched mean-field profiles are steeper(in 
units of 7 _1 ) than the ones observed, with better agreement as 7 is decreased. We believe 
that this discrepancy is due to the fact that T/T^J > T e ff/T®, as the equilibrium curve for 
7 > is flatter around T c than the mean-field curve. 

The need for a clear separation between boundary width and domain size may explain 
why in earlier molecular dynamics simulations a smaller exponent is found: in those 
computations the maximum domain size observed is only 6-8 times the range of the potential, 
so the true hydrodynamic regime was probably never reached. In simulations using the 
lattice Boltzmann method the interaction range is of the order of one lattice spacing and 
the hydrodynamic exponent is observed when the domain size is about 10 lattice units [f§], 
in agreement with our analysis. 



At early times the growth is consistent with a 1 3 behavior, for all potential ranges fl9| , |24 
This exponent is not associated with a scaling regime and it may not be universal, but it 
is not inconsistent with the experimental results (§,0. In fact, we were able to collapse the 
three curves onto each other through scaling of the lengths and times Fig. |J . 

Off-critical quenches were performed for a single range of the potential and three volume 
fractions. To our knowledge no published results of such simulations exist, although they 
are mentioned in J7|. Therefore, we compare our results to recent experimental work || 
and analyze them using the known coarsening mechanisms pP|J25] ] . In Fig. |3] we present the 



domain growth for 7 1 = 0.4 and volume fractions 0.16, 0.22 and 0.28. We plot R 3 (t).vs.t, 
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and the late times regime is clearly consistent with a | exponent. In this regime the system 
satisfies dynamical scaling very well. The presence of a late times ta growth at volume 
fractions less or equal than about 0.3 was observed recently in micro-gravity experiments 
|| and reasonably explained theoretically using a droplet coalescence mechanism p9| . |25| . 
This regime was also analyzed by Siggia who predicted a prefactor proportional to v^{v 
being the volume fraction of the minority species). In our simulations the prefactors are in 
reasonable agreement with the above prediction. Furthermore, we can clearly see the motion 
and coalescence of the droplets in movies of the dynamics. 

In writing down evolution equations for the order parameter and velocity fields in a sym- 
metric binary system with momentum and energy conservation it is generally assumed that 
the density and temperature variations are small |[26|| . We were able to check these quanti- 



ties directly by dividing the system into 'hydrodynamical' cells, each containing about 100 
particles, and computing the order parameter, density, temperature(fcsT = |m(( , u 2 ) — (v) 2 ), 
averages over the cell) and fluid velocity in each cell. While the statistical fluctuations of 
these quantities are not vanishingly small, as they should be for true hydrodynamic vari- 
ables, they are small enough for such a description to be at least reasonable. As already 
mentioned the total density is smaller at the boundary between domains. Away from the 
domain boundaries the values of the density and order parameter were close to their equilib- 
rium values. The temperature appeared to be uniform, with normal equilibrium fluctuations 
everywhere in the system. This shows that the time scales over which heat transport takes 
place are much smaller than the time scales over which there is significant phase separation. 

If, as argued, the linear growth regime is due to surface tension driven flows, one would 
expect bigger hydrodynamic velocities close to the interfaces than in the bulk. We looked 
therefore at the distribution of hydrodynamic velocities as a function of the order parameter 
for critical quenches. As the system segregates, bigger hydrodynamic velocities are typically 
observed close to the domain walls, identified by small values of the order parameter. While 
this may be due in part to the density being smaller at the boundary between domains, 
it is also consistent with the idea that at late times the flows are generated mainly by the 



curvature of well formed interfaces. We are planning to study this issue more carefully using 
bigger systems. 

A natural refinement of our model would be the use of the Enskog correction for the 
collisions [14j] (which would require numerically replacing the DSMC part with its recently 
proposed extension P7| , the Consistent Boltzmann Algorithm(CBA)) and of long range 
attractive forces between like particles. 

We believe that we have introduced a new model which is closer to reality than lattice 
gases usually simulated, but still tractable numerically by the new techniques that we have 
introduced. These techniques are of independent interest and have enabled us to compute 
temperature changes, density variations and profiles not done before. Furthermore, we can 
compare our results to some exact ones. In fact, we think that this model is a paradigm for 
the rigorous derivation of hydrodynamic equations. 

We would like to thank Frank Alexander for useful discussions. This research was sup- 
ported in part by AFOSR Grant 0159 and NSF Grant 92-13424. 
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FIGURES 

FIG. 1. Scaled two-point correlation function for critical quench with 7 _1 = 0.4. Final tem- 
perature T = 0.6T C °. 

FIG. 2. Scaled domain growth at critical quench; the 7" 1 = 0.3 and 7~ x = 0.5 curves have 
been collapsed onto the 7 _1 = 0.4 curve. A straight line fit(a = 1) is drawn for late times and a 
oc t3 fit (a = |) is drawn for early times. 

FIG. 3. Domain growth for off-critical quenches with 7"" 1 = 0.4; straight line fits are drawn. 

FIG. 4. Total density ra(*) and order parameter 4>(+) variations at the interface between the 
two phases; 7" 1 = 0.4 and T/T c ° = 0.6. The full lines are the solitonic solution(see text). 
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